ipo <- read.delim('ecol 563/ipomopsis.txt') ipo[1:8,] out1 <- lm(Fruit~Root+Grazing, data=ipo) coef(out1) # package for BUGS library(arm) # package for JAGS library(R2jags) # create variables x <- ipo$Root y <- ipo$Fruit # create dummy variable for grazing z <- as.numeric(ipo$Grazing)-1 n <- length(y) ##### BUGS/JAGS program stored in separate file: ipomodel.txt
model {
# likelihood
for (i in 1:n) {
y[i] ~ dnorm(y.hat[i], tau.y)
y.hat[i] <- b0 + b1*x[i] + b2*z[i]
}
# priors
b0 ~ dnorm(0, .000001)
b1 ~ dnorm(0, .000001)
b2 ~ dnorm(0, .000001)
# tau.y ~ dgamma(.001,.001)
tau.y <- pow(sigma.y,-2)
sigma.y ~ dunif(0,10000)
}
######## # set working directory where BUGS program file is stored setwd("C:/Users/jmweiss/Documents/ecol 563") # list of variables in model ipo.data <- list("n", "y", "z", "x") # initial value functions for chains ipo.inits <- function() {list(b0=rnorm(1), b1=rnorm(1), b2=rnorm(1), sigma.y=runif(1) )} # parameters whose posterior distributions should be returned ipo.parms <- c("b0", "b1","b2","sigma.y") # if WinBUGS program is in default location ipo.1 <- bugs(ipo.data, ipo.inits, ipo.parms, "ipomodel.txt", n.chains=3, n.iter=100, debug=T) # specify location of WinBUGS program ipo.1 <- bugs(ipo.data, ipo.inits, ipo.parms, "ipomodel.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=100, debug=T) # rerun with more iterations ipo.1 <- bugs(ipo.data, ipo.inits, ipo.parms, "ipomodel.txt", bugs.directory="C:/WinBUGS14", n.chains=3, n.iter=1000, debug=T) # JAGS run ipo.1.jags <- jags(ipo.data, ipo.inits, ipo.parms, "ipomodel.txt", n.chains=3, n.iter=1000) # WinBUGS components names(ipo.1) # JAGS components names(ipo.1.jags) # JAGS components in WinBUGS format names(ipo.1.jags$BUGSoutput) # model summary ipo.1$summary ipo.1.jags$BUGSoutput$summary # compare with frequentist results coef(out1) # samples of posterior distributions ipo.1$sims.matrix[1:10,] dim(ipo.1$sims.matrix)